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f— ^ , observed through partial and noisy data. In this paper, we consider 

l/^ • stochastic extensions in order to capture unknown influences (chang- 

0^ . ing behaviors, public interventions, seasonal effects etc). These models 

"^ ' assign diffusion processes to the time-varying parameters, and our in- 

CO , ferential procedure is based on a suitably adjusted adaptive particle 

C^ ■ MCMC algorithm. The performance of the proposed computational 

methods is validated on simulated data and the adopted model is ap- 
plied to the 2009 HlNl pandemic in England. In addition to estimat- 
ing the effective contact rate trajectories, the methodology is applied 
in real time to provide evidence in related public health decisions. Dif- 
; ^ \ fusion driven SEIR-type models with age structure are also introduced. 

Cu ■ population epidemic model; time-varying parameters; Bayesian infer- 

ence; Particle MCMC 

1 Introduction 

Epidemic models are often use d to simulate disease tr ansmission dynam- 



ics, detect emerging outbreaks (lUnkel and othersl . l2012l ). and assess public 



health interventions (|Boilv and othersl 120071 ). In order to capture the dy 



namics of epidemics, the main focus is generahy made on their intrinsically 
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dynamic elements such as the depletion of susceptibles or the population 
immunity evolution. Nevertheless, there are time-varying extrinsic factors 
that are crucial to the epidemic course. These may include social cycles 
(holidays), public interventions and climatic va riations. This has been il- 
lustrated for disea s es suc h as cholera, i nalaria (ICazelles and otherd . l2005l : 
lonides and otherd . l2006l ) or influenza ( Shaman and Kohnl . 120091 ). These 
studies were conducted eit her by relating climatic and incidence time-series 
([Cazelles and o^feersl . l2005l ) . which does not disentangle the effect of intrinsic 
and extrinsic factors, or by experimentally assessing the virus resistance in 
different climatic conditions ( Shaman and Kohnl . 120091 ) requiring an extrapo- 
lation to the population scale. Overall, the time-varying nature of epidemics 
poses a challengi ng statistical pr oblem stressing the need for suitable com- 
putational tools (JFergusonl . 120071 ). 

This paper considers a flexible modelling framework that encompasses 
time-varying aspects of the epidemic via stochastic differential equations. 
We aim at providing robust inferential procedures, incorporating the uncer- 
tainty associated with key parameters and accounting for data and model 
limitations. In order to provide an accurate and feasible computational 
toolbox, we provide Markov Chain Monte Carlo (MCMC) algorithms util- 
i sing recent developnients such as particle MCMC (PMCMC) algorithms 
( Andrieu and otfeer5l . l2010l ) and adaptive techniques ( Roberts and Rosenthal 
20091 ). Modelling aspects are presented in Section 2, while the computational 
framework is presented in Section 3. In Section 4 we evaluate the perfor- 
mance of the proposed adaptive PMCMC schemes on simulated data. In 
Section 5 we present various applications of the methodology to the 2009 
A/HINI pandemic, and conclude, in Section 6, with some relevant discus- 
sion. Further simulations can be found in the Supplementary Materials. 



2 Modelling framework 

2.1 Epidemic models with time- varying coefficients 

We adopt a SEIR model as a guide in this paper, although the methodology 
can be applied to other dynamical systems. The model is set in ([1]); S 
accounts for susceptible, E for infected but not infective, I for infective, and R 
for removed individuals. New infections occur at a rate /35f^, implying that 
the susceptible individuals make effective contacts at rate /3 (the effective 
contact rate), and only a fraction -^ of these contacts are made with infective 
individuals. The average period spent in compartments E and / is given by 
k~^ and 7""^ respectively. 
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The basic reproduction number, Rq, represents the number of secondary in- 
fections from a primary infected individual in a fully susceptible population. 
A related quantity is the effective reproduction number, Rt, refers to the 
number of secondary cases from an infected individual at time t. Rt is a, 
context-dependent quantity of high interest to policy makers as it indicates 
the possibility for the epid emic to grow {Rt > 1) or to decrease {Rt < 1) 
(JAnderson and Mavl . Il992l ). 

Epidemic models can be quite detailed (including individual character- 
istics, geographic information etc.) or basic, such as the SEIR model, that 
geographically aggregates the cases and assumes deterministic transmission 
processes, occurring at a given frequency each time infected and suscepti- 
ble meet. The latter are easier to estimate and interpret, but are based 
on strong assumptions that could lead to poor inference. In this paper we 
adopt stochastic extensions of the deterministic SEIR models. The addi- 
tional dynamic error is likely to contain structural mis-specifications and can 
subsequently be explored and potentially revised. We focus on large-scale 
epidemics, for which random effects in transmis sion processe s can be consid- 



ered to be well-approximated deterministically (JKurta . Il98ll ) . We adopt the 



paradigm that attributes the model limitations mainly to the time varying 
nature of the effective contact rate, henceforth denoted as /3t, rather than to 
the variability in individual characteristics or in transmi ssion processes. 



An early approach to estimate Rt can be found in iFine and Clarkson 
( 19821 ) ■ It can be implemented throug h discrete generation models or by re- 



const ructing the chain of transmission ( Cauchemez and o^/tersl . l2006l : lGriffin and others 



2011 



However, as Rt estimates contain both the effects of evolving trans- 
missibility and immunity, quantitative conclusions can hardly be generalised 
to situations where the immunological situation is different. We therefore 
concentrate on estimating /Jj rather than Rt- A number of approaches use 
a finite-dimension functi on space for the trajectory of /3f. Low-dimensional 
examples can be found in lCauchemez and other A ( 20081 ) . in which j3t is mod- 
el ed as a piece-wise linear functio n. In some higher-comp l exity models, as 
in lCauchemez and FergusonI ( 20081 ) and llonides and others ( 20061 ) . Pt is esti- 
mated freely with a few-weeks resolutions. Loosely speaking, as the number 
of parameters for the trajectory of Pt increases, model-induced biases fade 
out at the expense of the variance. A compromise is required to improve 
robustness and is often co i itrolle d through a regularising parameter. For 
example, in iHe and otherd (|201ll ). /3t is estimated using cubic splines, cali- 
brated via AIC 



2.2 Diffusion driven epidemic models 

We consider models where diffusion processes are used for some of the coef- 
ficients in ([1]). Although alternative formulations are possible, as discussed 



in Section I2.1| we focus on /3( to get 

\dxt = fJ.x{xt,0x)dt + a^{xt,6x)dBt, xt = h{fit), 

(2) 
where (J.x{^) denotes the drift, (Tx{-) the volatihty and h{) is a positive- valued 
function. The assigned diffusion may capture features such as behaviour 
changes, preventive measures, seasonal effects, holidays etc. When prior 
knowledge on /3( is available, it can be reflected in ^x{^) and (T^(-); e.g. if the 
contact rate is expected to converge, an Ornstein Uhlenbeck process can be 
chosen. Other options may include a sigmoid or a sinusoidal form; see for 



example (JRasmussen and otherm . l201ll ) . In absence of prior information or 



when the researcher wants to impose little restrictions, a Brownian motion 
can be used, with /Ux(-) = and crx{-) = cr (i.e. 6^ = cr). This model, with 
/i(-) = log(-), is henceforth denoted as BM. The obtained output can be ei- 
ther reported or used as an exploratory tool to construct a more structured 
model; see Section [231 for an application. The choice of BM implies a contin- 
uous, yet non-differentiable, path satisfying the Markov property. In cases 
where Pt is believed to evolve as a smooth function in time, higher order 
Brownian motions could be used. Loosely speaking, these may be regarded 



as eq uivalent to non-parametric approaches such as cubic splines (jWahba 



1990), with the model in ([2]) imposing a prior on /3( and a being a regu- 



larising factor. The rate /3f can be perceived as a product of a smooth and 
a rough component; the former being a population average of the intrinsic 
transmission procedure and latter containing extrinsic factors such as the 
amount of contact among individuals. It is therefore important to build a 
framework that contains both smooth and rough models. 

The above model can be estimated with an Extended Kalman Filter 



(EKF), as in lCazelles and Chad (jl997l ). EKF allows for fast computations. 



but is based on Taylor and Gaussian approximations whose error could be 
non-negligible; see Supplementary Materials for a relevant simulation ex- 
periment. Nevertheless, the EKF can still be used as a tool to construct 
efficient proposal distributions for MCMC schemes. It can also be used to 
optimize sequential Monte Carlo (SMC) algori thms, but either at a strong 



computational cost (ISarkka and Sottinenl . l2008l ) or crude time discretisations 
( Dukic and o^feersl . l2009l ) . Next, we develop a general framework for efficient 



MCMC schemes that allow for good approximations. 

3 Data augmentation via MCMC for diffusion driven 
epidemic models 

This section presents a general inferential framework for diffusion-driven epi- 
demic models. We adopt the Bayesian paradigm to incorporate parameter 



uncertainty and prior information in the estimates of /3t trajectories. The 
problem can also be cast as estimating par tially observed hyp o ellipt ic dif- 
fusions, thus presenting various difficulties ( Pokern and othera . |2009|). We 



begin by setting the model and justifying the need for data augmentation. 
Existing MCMC algorithms are considered but they can lead to extremely 
inefficient MCMC chains. We address the issue by taking advantage of the 
specific model structure to construct adaptive PMCMC schemes. 

3.1 Model and data augmentation setup 

For ease of exposition we focus on models satisfying ([2]) , but the framework 
covers models with different ODE systems or more time- varying coefficients, 
as in Section 15.31 Being in continuous time, t can take any value between to 
and tn- We denote the path of the ODE states vector Vt = {St,Et,It,Rt} 
between observation times ti and tj by Vi-j. The data, yi-n = {vti, ■■,ytn}j 
usually provide information for It at specific times (prevalence data) or for 
integrals of Vt (incidence data). In either case, we assume that they are 
obtained with error as the collection procedure is typically associated with 
additional uncertainty. The noise distribution is denoted with Pj^ with den- 
sity /(yi:n|Vo:n,^j/)- Note that, in the model of ([2|), Vt can be written as a 
deterministic function, g(-), of xt and the parameters 9^ = {k,-f,Vo). This 
function is the solution of the ODE and can be written as an intractable 
time integral involving xt- Hence, the model becomes 

\dxt = fix{xt,0x)dt + ax{xt,9x)dBt 

\yi:n\Vo:n, Oy ~ Py(yi;n| Vb;„, 9y), Vq-u = g{xO:n,9v) 

Denote with F^ the distribution of the diffusion xt defined from the SDE 
above. We require the existence of a unique weak solution which translates 
into some mild assumptions on ^x{-) and a.rf.); e .g. lo cally Lipschitz with a 



linear growth bound, see for example lOksendall (J2003|). The distribution of 



¥x may also be viewed as a prior on xt, or else Pf The model can now be 
defined from ¥y, P^, and the assigned priors on 9 = {9y, 9v,9x}, denoted by 

7t{9) 

T^{x0:n,9\yi:n) OC /(yi:„|Vb:„, 9y) X dP^. X Tt{9) (4) 

Given direct observations on xt , it would ha ve been possible to draw ap proximation- 



free inference on dFx using the approach of lBeskos and otherd (|2006l ). How- 
ever, this is not possible in our case given the non-linear functionals in g(-) 
that render (3.4) intractable. We proceed by discretizing the path of xt, and 
therefore of f3t and Vt- More specifically, we introduce m points between 
each pair of successive observation times ti and tj+i (i = 0, 1, ... ,n — 1). 
When referring to the discrete representation of a path, the superscript dis 
will be used; for example for a step S = ^^^jry, the discrete skeleton of xt will 



be denoted by Xq.^ = {xq, X5, X25, • • • , 3^t,i}- The presence of x^.f^ allows for 
approximations of ^ through the Euler-Maruyama scheme to evaluate dFx 



iS<tnP(^is\X(i-l)S,0x), 






Moreover, given Xg*^,the ODE can be solved numerically to obtain Vq^^'* and 
evaluate /(•). The approximation error can be made arbitrarily small by 
increasing the user-specified parameter m. 



3.2 Data augmentation via Gibbs schemes 

Mode l (l3l) can be put in th e conte xt of lChib and o^/teral ( 20061 ) . iGolightlv and Wilkinson 



(|2008h 



or 



KalogeropoulosI ( 20071 ). In these approaches, a Gibbs scheme can 

The data 



be used to sample from the joint posterior in (JU of x^*^ and 6. 



augmentation algorithm alternates between drawing Xq*^ given 9, and updat- 
ing 9 conditional on the augmented path Xq*^. The MCMC protocol ensures 
that the chain provides samples from the marginal posteriors of Xq*^ and 
9. Nevertheless, the properties of the algorithm may become unacceptably 
poor. There are two essential issues associated with such schemes. The first 
concerns the non-trivial step of sampling on the diffusion pathspace of Xf. 
The second problem is caused by the high posterio r correlations between 

x: 



dis 



and 0, leading to reducible chains as m increases ([Roberts and Stramer 



2001 



The majority of the literature on data augmentation schemes for diffu- 
sions handles the conditional updates of Xq'* with an independence sampler. 
As it is difficult to find good proposal distributions for the entire Xq*^, the 
path is usually split into blocks. Overlapping blocking strategies are essential 
to ensure that all points are updated and continuity of the path is retained. 
An alternative way to update Xn^^ is to u se the particle filter via the Particle 
Gibbs algorithm of lAndrieu and otherd ( 2010l ). But unless the issue of high 



posterior correlation between Xq*^ and 9 is resolved, none of these schemes 
will improve the overall MCMC performance. The problem is caused by the 
quadratic variation process of xt that identifies 9x- For crx{xs,9x) = a we 
get 



lim 

5^0 



(6) 



y^ (xis - X(^i_i)s)'^ = / a'^ds = a'^{tn - to) 

to<i5<tn ° 

Thus, the conditional posterior of a converges to a point mass as 5 tends to 0. 
In practice this translates int o an increasingly slow MCMC algorithm with 
a conver gence rate of 0(rn) (IRoberts and Stramen . I2OOII ). Schemes with a 
fixed m ( Cori and ot/iersl . l2009l ) could work in some occasions but the approx- 
imation error could be substantial. In some cases, the problem can be tack- 
led with suitable reparametrisation. The approach of [Roberts and Stramer 



(|200ll ) involves transfor ming to a diffusion xt wi th unit volatility. An alterna- 
tive scheme is offered bv lChib and other^ i 20061 ) where the driving Brownian 



motion of xt is being used. In these algorithms the ODE states vector Vi^!^ 
becomes a function of a, XQ-n and 9^. Hence, in a Metropolis step, every 
proposed value of a* is associated with the corresponding values of Vij:^^ . 
This succeeds into breaking the perfect dependence between Vq:^^ and a, 
even for m — )■ cx). But since components of V^l^ (or functionals thereof) are 
observed with error, the associated proposed values Vf^l^ should be close 
to the data for the move to be accepted. As the observation error becomes 
small and the data increase, this becomes increasingly difficult and leads to 
very small moves for a and poor MCMC mixing. More details and simula- 
tions supporting this argument are provided in the Supplementary Materials 
(Appendix E). Consequently, we overcome this issue by updating x^:^ and 
9 jointly via the PMCMC algorithm, which is essential as it is not straight- 
forward to implement joint updates with the other approaches mentioned in 
this section. 

3.3 Adaptive Particle Markov Chain Monte Carlo algorithms 

Particle filters are SMC algorit hms used to recursive l y exp lore conditional 



densities in state space models (jDoucet and Johansenl . 120111 ) . For given val 



ues of 9, N particles (x^) are sequentially propagated from to to tn- In var- 
ious time steps tj, the trajectories that best fit the data yi-i are given more 
weight through resampling. Algorithm [1] shows how they can be applied in 
our context. The quantity L*"*"^(^) provides unbiased estimates of p{yi-i\9) 

Algorithm 1 Particle Filter algorithm 

Initialise: Set L'^{9) = 1, W^ = jj^, sample {xq)j=i,...,n from p{xo\9) and 
calculate (V^"')j=i,...,Ar by solving the ODE (for example with the Euler 
scheme) 

for i = to n — 1 do 
for j = 1 to A^ do 

Sample (x^.j , ]^) from ([5]) and calculate {V^.^,-^^) by solving the ODE 

Set a^ = f{y^+i\Vl^,) 
end for 

Set W^, = ^,^, and U+^O) = L\9) x i E «^' 

Resample (Fo^.j+i,x;J.j+i)j=i,...,Ar according to (W/+i), 
end for 

and the resampling step is essential to control the variance of that estimate 
over time. Algorithm [1] also provides a random sample from p{xi-i\yi:n,0)- 
In order to sample from 'iT( x-i;n,9\y-i;„), the PMCMC a lgorithm can be used. 



PMCMC was introduced in lAndrieu and other ^ (2010|) and successfully inte^ 



grates particle filters in MCMC algorithms. Its implementation is presented 
in Algorithm [2j The issues of Section 13.21 are now addressed as xjj*^ and 9 

Algorithm 2 Particle MCMC algorithm (particle Marginal Metropolis Hast- 
ings version) 

Initialise: Set current value, 0, to an initial value. Use Particle 

Smoother (PS) according to Algorithm [1] to compute p{yi:n\G) = L[9) and 

sample x\.^ from p{xi;n\yi:n-, 0) 

for It = 1 to N Iterations do 
Sample 6* from Q{e, .) 

Use PS to compute L{9*) and sample x\.^ from p{xi:n\yi:n)0*) 
Do ~e = 9* (and xf^„ = xfj with probability 1 A ^^^g^gg-f^^ 

Record 6 and x^.^ 
end for 

are sampled jointly. In other words Xg*^ is being numerically integrated out, 
while a sample from its posterior is obtained at each MCMC iteration. 

While the PMCMC algorithm is theoretically valid even for a single par- 
ticle, large values of N are usually required for reasonably stable accep- 
tance rates and large moves in the 9 space; see the Supplementary Mate- 
rials for a relevant simulation exercise. It is therefore essential to update 
the d-dimensional 9 at once, marking the proposal Q{9, .) crucial to the 
overall MCMC performance. In this paper we propose to use the adap- 
tive Metropolis algorithm of [Roberts and Rosenthall ( 20091 ) . After trans- 



forming the parameters to take values in the real line we use a Normal 
distribution centered at the current value of 9 and with covariance given 
by eE. Static random walk metropolis proposals set T, = I^ or T, = 'S 
and tune e to obtain acceptance rate of 0.234. Adaptive schemes change 
the value e for each iteration i through diminishing adaptation; e.g. by 
ej+i = exp {log(ei) + a"(AccRate - 0.234)} where ai = 0.999 and 'AccRate' 
denotes the acceptance rate up to iteration i. The covariance matrix Sj+i 
can also be updated as 

a2M (^9, e^So) + (1 - a2)M (^, 6^S, 

where a2 is usually set to 0.05, Sj is the posterior covariance matrix es- 
timated by the draws up to i and Sq should be specified in advance. In 
this paper we enhance the above adaptive algorithms utilising information 
from the EKF to estimate the covariance S or T,q. One choice, EK-Mode, 
is the observed information matrix at the mode identified by EKF, evalu- 
ated through numerical differentiation. Another choice, EK-MCMC, is to 
run an approximate MCMC scheme based on the EKF approximation of 
the likelihood and compute the posterior covariance from the draws. Note 



that the computational burden of these methods is marginal with regards to 
the PMCMC. As demonstrated in Section [H the use of EKF can result in 
substantial improvement. 

4 Simulation Experiments 

The proposed algorithms are illustrated and tested on simulated data in 
this section. We focus on the BM model, where log{f3t) follows a Brownian 
motion with volatility a, corresponding to the case of little information on 
the shape of (3t- The trajectories of (3t were drawn either from the BM model 
itself (experiment 1) or from a deterministic sigmoid curve (experiment 2). 
The data yj, i = 1, . . . , 50 represent noisy observations of weekly new cases 
of the epidemic J^^^p. ^ kEtdt. We complete the model by assigning a Normal 
distribution to each log(yj) with mean \og{ j_^^^^^kEtdt) and variance r^. 
The parameters were tuned to obtain realistic epidemic incidence curves, 
and observations were generated setting r = 0.1. The assigned priors were 
informative for A;, 7 and -R(to) ^-nd vague for E{tQ), /(to), cr and r, as in 
Section O We used 3,000 particles and 100,000 MCMC iterations after 
a long burn-in period. Fig. [1] shows estimates and 95% pointwise credible 
intervals of the path, provided by the adaptive PMCMC initialized with 
EK-MCMC. The posterior output is in good agreement with the simulation 
trajectories suggesting that the underlying trajectory of Pt can be estimated 
reasonably well from the partial and noisy observations considered. More 
can be found in the Supplementary Materials (appendix C), where we also 
considered a value of r = 0.05 and obtained similar results. 

Next, we use the data of experiment 1 to compare the proposed adap- 
tive PMCMC schemes. Comparison is made in terms of the effective sample 
size ESS = (1 -|- 2^^>-^ 'n(i))~ \ with J^,^ r]( i) being the sum of the lagged 



sample auto-correlations, as in iGeveiJ ()1992l ). We record the minimum ESS 



among the MCMC components and multiply by 100 to monitor the per- 
centage of the total iterations that can be considered as independent. We 
consider three covariance matrices for each of the two adaptive algorithms 
defined in Section [3.31 I^ and the ones from EK-Mode and EK-MCMC. For 
the schemes that adapt e the minimum ESS was 0.008% (Id), 0.19% (EK- 
Mode) and 0.54%(EK-MCMC), whereas for the schemes that adapt S we got 
0.57%, 1.24% and 1.38% respectively. Clearly, adapting E is crucial to obtain 
a reasonable performance, unless the matrices obtained from EK-Mode or 
EK-MCMC are used. The proposed adaptive algorithms induce substantial 
improvement that is expected to intensify as the dimension of 9 increases. 



5 The 2009 A/HINI pandemic 

5.1 Data, model and estimates 

The proposed methodology is iUustrated on data from the A/H1N1(2009) 
pandemic in England between June and December 2009. The data consists 
of estimates of weekly ILI c ases yi-n given by the Health Protection Agency 
( Baguelin and o^/iersl . l2010l ). The estimates were obtained from the recorded 



ILI cases among a selected sample of GPs. They accounted for over-reporting 
due to similarities in symptoms with other respiratory diseases, based on sub- 
sequent virological positivity tests. Corrections for asymptomatic infections 
and the patients' propensity to consult were also made. Overall the two 
datasets are different by a multiplicative coeffi cient c = 10, whose val ue is 
also supported by a further serological survey ( Miller and other A . |2010|). In 



our analysis c is initially held fixed to 10, but this choice is explored fur- 
ther in Section 15.21 We adopt a model that admits noisy data to reflect the 
associated uncertainty. The noise model of Section 2] was used, combined 
with a BM formulation of P^;. Vague priors, A^>o(0, 10^), were put on r, a 
and /3n- The priors for k an d 7 were obtained from additional data sources 
( Baguelin and otherd . |2010| ). the results of which are summarised through 



Normal distributions that place 95% probability in a symmetric manner be- 
tween 1.55 and 1.63 days for the latent period k~^, and between 0.93 and 
1.23 days for the infectious period 7^^. A Dirichlet distribution was used for 
the initial proportions in compartments S, E, I, R, constraining the mean of 
the one in R to be 0.15, its variance 0.15^, and the means of the other initial 
proportions to be equal. 

The adaptive EK-MCMC algorithm was applied to the data and Fig. [2] 
depicts the incidence curve together with the posterior mean and pointwise 
95% credible intervals. Estimates of /3t are also displayed indicating various 
changes over time. The changes in Pt are consistent with the argument that 
schools closure for holidays have been driving the epidemic: different values 
are observed during school and holidays periods, appearing to be synchro- 
nised with schools opening and closing. Posterior summaries for the static 
parameters, as well as a sensitivity analysis on the priors can be found in 
the Supplementary Materials. These suggest that inference is quite sensitive 
to the choice of prior for k and 7, but not for the remaining parameters. 
It would be interesting to repeat the procedure under an evidence synthesis 
framework and vague priors. 

5.2 Application in real time. Was the first wave waning due 
to depletion of susceptibles? 

In this section the methodology of the paper is applied in real time, i.e. 
considering partial datasets from June 2009 up to the 20th of July, the 7th 
of September and the 26th of October. Each time the algorithm is run from 

10 



scratch to provide samples from the joint posterior iT{xi-i,9\yi-i). From a 
computational cost point of view this procedure can be improved further 
by utilising previous MCM C runs, for example under the SMC^ framework 



(jChopin and otfeersl . l201ll ). We did not pursue this direction further, as the 
PMCMC algorithm runs quite fast (less than 2 hours on a standard PC). 
In order to reduce uncertainty, especially at early stages, the value of r was 
set to 0.1 rather than being estimated as in Section [5.11 We otherwise use 
the same model as before. A model with integrated Brownian motion was 
also fit but BM was chosen in terms of DIG; see Supplementary Materials 
(Appendix C). The main results are shown in Fig. [3l 

On August 1***, the first wave of the epidemic had waned, incidence rates 
were decreasing and schools had closed. There were two competing scenarios 
to explain the epidemic decline: (i) holidays had caused the waning of the 
epidemic by lowering the effective contact rate. Hence, a similar or stronger 
wave could occur when schools would reopen in September in colder climatic 
conditions, (ii) The epidemic had stopped independently of holidays because 
a critical proportion of the population had been infected, conferring a suf- 
ficient level of herd immunity to stop the epidemic. In this case, no second 
wave was to be expected in Se ptember. On August 1^* th ere was great un- 
certainty around the value of c ( Baguelin and o^/iersl . 120101 ) . which is crucial 



in distinguishing between the two scenarios. We therefore conducted the 
following exercise. 

The PMCMC algorithm, run up to August 1st, provides samples from 
the posterior of the difference in /St between July IS*'^ (before the decrease 
in incidence) and August 1***. For c = 10, the 97.5% point of this posterior 
is —0.32, indicating a decrease in /3j. The latter supports scenario (i), as the 
competing scenario is associated with a zero-decrease in (3t- Nevertheless, 
as this value depends on c, the algorithm was run for different values of it 
ranging from 20 to 150. The results appear on Fig. 31 Note that the 97.5% 
point of interest increases as a function of c and reac hes for a correction 



facto r close to 70. As this level seemed unrealistic ([Baguelin and others 



2010l ) , the experiment provides evidence in favour of scenario (i) highlighting 



the danger of a second wave in September, that actually occurred. Such 
evidence can be important for decision-makers, especially when considering 
implementations of preventive measures as vaccines. 

5.3 A multiple age group diffusion driven SEIR model 

The analysis of Section [5T] can be used to construct more structured models. 
For example, the effect of holidays is evident and may differ from children to 
adults, thus casting doubts on the assumption of a homogeneous population. 
It seems more natural to consider a model with two age groups (c:children 
and a:adults) and target all possible effective contact rates among them. In 
our notation f3'^"' refers to the effective contact rate from children to adults 
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and S"^ denotes the number of susceptible children. For reasons of parsimony 
we assign Brownian motions to log{/3f'^), log(/3f") and treat /3'^, 13'^^ as con- 
stant. We also set we set l]^"" = 13"''^ = b, i n line various multiple ag e grou ps 



epidemic models in different settings (e.g. IWhitaker and FarringtonI (J2004l ))- 



The dynamic part of the model is now given by 

cttb^ Q(2 I Qcc ^t \ P\ t \ ""^t Qc { pec ^t I Ix -^i \ h, TpC 

~'-'t \Pt TI^ + "J!^ I ) ~dr ~ '-'t \Pt Arc "•" Ojttj: I — Ki/j , 



dt ~ '-'t \^t Arc ^ "N<^ J ) dt ~ '-'t yyt Arc t ^ j^a 



dt ~ '-'t yi^t N'^ ' ''N'= J ' dt ~ '-'t yi^t N'^ ' ""N' 



The data from the A/H1N1(2009) pandemic provide incidence estimates for 
children and adults separately so they can be used to estimate the model 
of d?]). If only final outcome data were available, not all effective contact 
rate parameters would be estimable. However, the temporal dataset pro- 
vides extra information by the relative variation of susceptible and infective 
population in adults versus children. We applied the EK-MCMC scheme, 
which was essential in order to obtain reasonable MCMC performance. Fig. 
[5] depicts the results. Unlike earlier attempts with versions of a multi-group 
model with a single diffusion driving all contact rates, the fit appears to be 
good. The trajectory of children seems to be similar with that of Fig. [2] thus 
stressing their role to the evolution of the epidemic. More details, including 
posterior summaries for the parameters and information about the priors 
can be found in the Supplementary Materials (Appendix C). 

6 Discussion 

In this paper we examined epidemic models where some of the parameters are 
represented by diffusions or integrals thereof. The main motivation was to 
account for various time varying drivers (virus evolution, seasonality, schools 
closure, etc), while maintaining a simple interpretation. We present a unified 
framework that supports data augmentation MCMC schemes based on fine 
partitions on the diffusion path. The associated approximation error can be 
controlled by the user without affecting the M CMC performance and can b e 



viewed as an extension of the approaches by [Roberts and Stramen (|200ll ) ; 



Chib and otherA ( 20061 ) to the more challenging observation regime of this 



paper. The consideration of the algorithms in a continuous time setting re- 
vealed major issues associated with Gibbs data-augmentation schemes. This 
justifies the use of particle MCMC, which updates paths and parameters 
jointly, while pointing directions for future research on Gibbs schemes. We 
also presented a cornputati onal machinery based on the PMCMC algorithm 
( Andrieu and o^/tersl . 120101 ) . that was integrated in an adaptive MCMC con- 



text. We consider EKF based adaptive algorithms that can offer substantial 
improvement, especially in cases with many static parameters. This paper 
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is one of the first appHcations of PMCMC in epidemic models and data ; 



standard PMCMC schemes were also used in iRasmussen and otherd (|201ll ). 
Initially we relied on a simple SEIR model but such an analysis can be 
viewed as an exploratory tool towards more structured models; e.g. the 
age-structured model of Section 5.3 that appears to be an improved rep- 
resentation of reality. This approach can help in developing richer mod- 
els and testing alternative scenarios for public health interventions, or to 
bring further insights on extrinsic factors such as climate on the dynam- 
ics of epidemics. Moreover, this frame work can support multiple so urces 



of data, of potentially different nature: IRasmussen and otherd (|201ll ) has 
shown how time series and genealogies can be combined in a PMCMC in- 
ference framework for more informative estimates. While we worked mainly 
with infiuenza time series, the developed methodology can be applied to other 
cases; current work cons iders its application as part of the CHARME project 



(JBoilv and otherd . 120071 ) . The presented approach may also be thought as 
an alternativ e to the white no i se mo deling of environmental stochasticity 
introduced in lBreto and otherA ( 20091 ) . as it offers to the possibility to cap- 



ture the dynamics of environmental drivers. A potential next step will be to 
combine environmental with demographic stochasticity, modelling infections 
as Poisson processes which rates depend on a time-varying j3. 

The inferential framework presented in this article shares the "plug and 
play" feature of the Iterated Filtering methodology. While extra care and fur- 
ther study is required for specific models or datasets, its algorithmic aspects 
can be decoupled from the modeling aspects. This provides the possibil- 
ity to develop generic inference packages: we are currently working towards 
its integration in a generic inference platform inspired from the R package 
POMP. 

7 Supplementary Materials 



Supplementary material is available online at http://biostatistics.oxfordjournals.org 
It contains implementation details for the PMCMC algorithm (Appendix A), 
a comparison with the EKF (Appendix B), additional information on Sec- 
tions [Hand [5] (Appendix C). Also the sensitivity analysis (Appendix D) and 
a detailed exposition of the issues of Section [3.21 (Appendix E). 
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Figure 1: Illustration of how the underlying dynamic of the effective contact 

rate can be estimated from weekly recorded cases. 

Green dots indicate simulated observed incidence (top panels). Green lines 

indicate simulated effective contact rate trajectories (bottom panels). Black 

dotted lines indicate the mean of the pointwise posterior density. Dark and 

light blue areas show credible intervals, respectively at 50% and 95% levels. 

Top panels: simulated weekly numbers of cases observed with noise, and 

corresponding model-based offline reconstructions (left: experiment 1, right: 

experiment 2) 

Bottom panels: simulated and estimated trajectory of the effective contact 

rate (left: experiment 1, right: experiment 2) 
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Figure 2: Weekly incidence data from the A/HINI 2009 influenza pandemic 

and corresponding offline estimates of the effective contact rate. 

Green dots indicate incidence estimates provided by the Health Protection 

Agency. Black dotted lines indicate the mean of the pointwise posterior 

density. Dark and light blue areas show credible intervals, respectively at 

50% and 95% levels. Holidays are indicated by a light grey area. 

Top: observations of the weekly total number of A/HINI influenza cases in 

London (per 100 000 inhabs.) and model-based offline reconstruction 

Bottom: offline estimates of the effective contact rate. 
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Figure 3: What could have been inferred by carefuhy fohowing the epidemic 

in real time? 

Green dots indicate observed incidence estimates provided by the Health 

Protection Agency (left panels). Black dotted lines indicate the mean of 

the pointwise posterior density. Dark and light blue areas respectively 

indicate 50% and 95% credible intervals of the posterior density. Holidays 

are indicated by a light grey area. 

Left panels: HPA estimates of the weekly total number of A/HINI influenza 

cases in London (per 100 000 inhabs.) 

Right panels: "real-time" estimates of the effective contact rate. 
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Multiplicative correction coefficient c 



Figure 4: The implication of different scenarios for the real value of under- 
reporting on the decrease of the effective contact rate between July 13* and 
August I'** 

For each value of c, the mean of the posterior density for (3(August f**) — 
(3{July 13*^) is plotted in black. Dark and light blue areas respectively in- 
dicate 50% and 95% credible intervals of the posterior density. The dotted 
line locates the scenario whith no change in the effective contact rate. 
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Figure 5: Offline estimates of tlie effective contact rate among cliildren 

and adults during the A/HINI 2009 influenza pandemic using a 2-classes 

age-structured model and age-specific incidence data. 

Green dots indicate observed incidence estimates among each age group 

provided by the Health Protection Agency (first and second panels). Black 

dotted lines indicate the mean of the pointwise posterior density. Dark and 

light blue areas respectively indicate 50% and 95% credible intervals of the 

posterior density. Holidays are indicated by a light grey area. 

First panel: HPA estimates of the weekly total number of A/HINI influenza 

cases among children in London (per 100 000 inhabs.) 

Second panel: HPA estimates of the weekly total number of A/HINI 

influenza cases among adults in London (per 100 000 inhabs.) 

Third panel: offline estimates of the effective contact rate from children to 

children. 

Fourth panel: offline estimates of the effective contact rate from adults to 

adults. 
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Appendix A: details of the PMCMC implementation 

In this Appendix, we provide more details for the practical implementation 
of the PMCMC algorithm presented in this article. We specify how to deter- 
mine key parameters of the algorithm, i.e. the Euler discretisation time-step 
6, the number of particles Nparts used in the Particle Smoother (PS), and 
how to set the Metropolis updates of the parameter vector 9. 
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Determining the Euler discretization time-step 

In general, solutions of the nonlinear ODEs encountered in epidemic models 
are not available in closed form. In order to evaluate iT{xQ:n,0\yi-n), the tra- 
jectory of the system needs to be discretised according to a given time-step 5 
to provide an approximate solution. The Euler approximation ensures that 
as 5 tends to 0, TTs{xQ-n,0\yi:n) converges to '7r(xo:n) ^|2/i:n)- In practice a se- 
quence of decreasing 5 values is chosen and quantities such as E[ps{a\yi-n)] 
or E[ps{T\yi-n)], are monitored. For sufficiently small values of 6, a conver- 
gence is generally observed, as shown in Fig. [T]for two different datasets of 
weekly influenza data. In this case, a, 6 = 0.1 day would be a reasonable 
choice. We note at this point that the computational cost is of 0{5^^) 

Determining an optimal number of particles 

The PMCMC algorithm is theoretically valid regardless of the number of par- 



tidesNparts used in the particle smoother, as shown in lAndrieu and others 



( 20101 ). Nevertheless, the smaller the number of particles used in the particle 



smoother, the more noisy the estimate of the likelihood p{yi:n\(^) becomes. 
This noise has a negative impact on the acceptance rate of the MCMC al- 
gorithm run in the 9 space. Consequently, Nparts has to be big enough so 
that it won't affect the acceptance rate, keeping in mind that the cost is of 
0{Nparts)- Hence, a compromise needs to be achieved. Fig. [2] shows how 
the acceptance rate increases as the number of particles gets higher. Note 
that the acceptance rate has a plateau form, indicating that it is perhaps 
not worthwhile to increase Nparts beyond some point. We repeat the exper- 
iment for two different values of the measurement error parameter (r) . This 
figure shows that when the observational noise decreases, more particles are 
needed, which is explained by the fact that the particles need to be 'closer' 
to the observations. 

Appendix B: Assessing the validity and limitations 
of the Extended Kalman approximation 

Here, we compare the particle filter with the EKF approach to assess the the 
extent of the gain of avoiding the approximations of the latter. A set of 100, 
7-month long, time-series of weekly influenza cases were drawn from the BM 
model. In order to ensure realistic epidemic datasets, we 'reverse-engineered' 
randomly selected influenza time-series (y].'^^''')^^^ fro m the freely avail- 
able Google FluTrend data ( Ginsberg and others . l2008l ). For each of the 



datasets, we obtained estimates of (/3q.^'"' ) ^2!i and the corresponding pa- 



rameters {Oj)j^^. These quantities were then used to generate influenza 

)" 



time-series (y^.^'-')y2?i- The static parameters of the model (initial condi- 
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tions, k, 7, a and r) were assumed known, to isolate the problem of es- 
timating /3i from accounting for parameter uncertainty and perform more 
relevant comparisons. We compare the following two estimators of /3^ *'"'"' : 
;5j * '"' = £'(/3^|y|.j, 0p obtained from the filtering distribution where E{.) 
denotes Monte Carlo estimates of the relevant expectations, and the the 
EKF estimator /3f^^'J' = E^^^{j3l\y{.-,0*) where E^^^{.) denotes expec- 
tation under EKF. The performance of the estimators is measured through 
their bias and Mean Squared Error (MSE). The results indicate a better 
performance for /3j * '■'. The bias of the estimates provided by the EKF is 
0.0285 while use of /3f^'*'-'' reduces the bias by about 78% (0.0063). The 
corresponding relative reduction in MSE is smaller (10%, 0.0270 to 0.0242), 
indicating a bias- variance tradeoff. Use of the smoothing distribution esti- 
mator /3f™'^' = E{l3J\y^^^'^,e*) is associated with a further 87% (0.0032) 
reduction in the MSE, while keeping the bias at the same low levels. The 
estimators f3- * '-' and /3- '"'"' are associated with a tolerable computational 
cost of 2 hours on a standard PC. In conclusion, the bias introduced by the 
Extended Kalman approximation is non-negligible with regards to the level 
of accuracy that can be obtained with exact particle methods on this type 
of datasets. Nevertheless, this study has shown the approximation to be 
robust and motivates the use of the approximated model as a proxy for the 
exact one, for example to initialize the particle MCMC algorithm in the way 
presented in this article. 



Appendix C: details on the simulations and results 
of sections 4 and 5 

Assessing the validity of the particle MCMC through simula- 
tions 

This section deals with a series of experiments that were conducted in order 
to assess the validity of the particle MCMC, and to illustrate on different 
examples how the trajectory /3i:„ could be captured from noisy weekly cases 
observations. Experiment 1 is based on a trajectory of /3i;„ simulated from a 
random walk model with volatility cj^ = 0.07^. Two different corresponding 
epidemic datasets have been generated, for given and equal initial conditions 
and biological parameters, respectively with observational noise r = 0.1 
(experiment l.a) and r = 0.05 (experiment l.b). Similarly, two epidemic 
datasets resulting from an effective contact rate following a significantly de- 
creasing sigmoid were generated with respectively r = 0.1 (experiment 2. a) 
and r = 0.05 (experiment 2.b). 

For each of these datasets, our proposed methodology was run to estimate 
a, T and /3i:„. cFig. 1 of the main text contains estimates of the incidence 
time series as well as f3i-n for the experiments with r = 0.1 (Exps l.a and 



Capturing the time-varying drivers of an epidemic 4 

2. a). Corresponding figures for Experiments l.b and 2.b are shown in Fig. 
[3] of this document. Moreover, Table [2] presents the mean, median and 95% 
credible intervals for the estimates of a and r in each of the experiments. 
In experiment 2 the estimates seem to be in good agreement with the true 
values, as the latter are contained in the 95% credible intervals. The aim 
was to assess the robustness of the proposed methodology to model mis- 
specification, fitting a Brownian motion to a smooth sigmoid curve. The 
algorithm performs reasonably well, succeeding in capturing the trajectories 
of f3t and all the parameters except for r which is slightly underestimated. A 
potential explanation for this is that part of the variability is absorbed from 
the volatility parameter of the Brownian motion. 

A/HINI pandemic 

We provide here with the corresponding trace plots (Fig. [3]) for the param- 
eter estimates of Section 5.2. The trace plots indicate good mixing of the 
PMCMC algorithm for every parameter. This was achieved by following the 
procedure presented in Appendix A. 

Age-structured model: children and adults 

Additional information on the analysis of section 5.3. can be found in Table 

El 

Illustration of alternative approaches on the real time example 

In Fig. [5l we repeat the real-time analysis, conducted in the Section 5.2 
of the main paper, under two alternative approaches. First, we consider a 
model with an integrated Brownian motion (IBM) on xj, implying smoother 
/3t trajectories as opposed to the non-differentiable paths induced by the 
Brownian motion (BM) formulation. The choice between those models is 
not trivial and could depend on the context of the epidemic. We decided to 
adopt the model with Brownian motion (BM) on the bas i s of t he Deviance 
Information Criterion (DIG) of ISpiegelhalter and others] ( 20021 ): as we can 



see from Fig. 3 of the main document and Fig. [5] of the SM, the BM model 
is consistently better in that respect. Nevertheless, there is uncertainty on 
the performance of the DIC criterion in the setting of this. It would be quite 
interesting to explore this further in the future and compare with alternative 
model choice criteria. 

The second approach adopts and applies the metho dology of maximum 



likelihood via iterated filtering (MIF), introduced in (jlonides and others 



20061 ) . Initially, estimates 9 are obtained by maximising p{yi;n\&) subject 



to some constraints set by the priors. Second, a particle filter was run with 6 
fixed to its estimated value. As expected, given that we are not accounting 
for parameter uncertainty, the resulting pointwise 95% credible intervals are 
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narrower; roughly 50% on the 6-month dataset and even more at early stages 
with less information on 6. 

Appendix D: Sensitivity analysis 

We explore in this section the robustness of the obtained estimates. We con- 
centrate on the example of Section 5.1 in the main paper, where the obser- 
vational noise was estimated and informative priors were used for the initial 
proportion of immune individuals in the population (i?(0)), and the lengths 
of the latent {k~^) and infectious periods (7^"*^ ). For each of these quanti- 
ties, the mean of the prior densities have been tilted by respectively -20%, 
-10%, +10% and +20%. In all cases but for the tilted parameters them- 
selves, the median of the corresponding posterior densities lie both in the 
95% and 50% credible intervals of the untitled case. The resulting medians 
estimates are shown in tabled along with the originally obtained summary 
statistics. Furthermore, for all the runs with tilted priors, the resulting /3i-n 
trajectories remain completely in the original 50% credible intervals. 

Appendix E: details of Gibbs data augmentation scheme 

In this section we provide more details on the Gibbs schemes discussed in 
section 3.3 of the main text. Stochastic epidemic models presented in this 
paper can be written as 

\dxt = fixixt,0x)dt + ax{xt,9x)dBt, < t < tn ,. 

where Vt represents the ODE states vector observed trough partial and noisy 
data yi:n- The rest of the model is defined in section 3.1. Since it contains 
intractable densities we work with the time discretised versions j;g** and 
V^l^ and proceed using the Euler approximating scheme. A Gibbs algorithm 
alternates between updating the trajectories of Xq*^ (and consequently V^l^) 
given 9, and vice versa. Nevertheless, as the Euler time step 5 goes to 0, the 
quadratic variation process of Xt uniquely determines the va lue of 9x in ax{-) 



and the algorithm degenerates (jRoberts and Strameij . l200ll ) . In practice this 



translates into a mixing time of 0{m). 

In the context of diffusion driven epidemic models t his problem was 



dealt wit h suitable reparametr isations such as the ones in IChib and others 



( 20061 ) or lKalogeropoulod ( 20071 ) . The latter uses the Lamperti transform, i.e. 



Xt -^ H{xt,9x) = ri{xt,9x) — 'ri{xQ,9x) ='■ Ut where r]{-;9x) is an antideriva- 
tive of a~^{-;9x)- Assuming that ax{-',9x) is continuously differentiable, an 
application of Ito's lemma provides the SDE of the transformed diffusion ut 
as: 

dut = y{ut] 9x)dt + dBt , uq = , (2) 



Capturing the time-varying drivers of an epidemic 6 

where 

Let P** denote the distribution of ut- Girsanov formula provide its density 
with respect to that of a standard Brownian motion, denoted by W, 

,m>u f rtn -^ rtn ^ 

expM v{us]6x)dus- - I v{us]exfds\, (3) 



dm 

and the state vector can be written as 

VQ:n = K{uQ:n,Xo,9x,9^). (4) 

The model can be defined from Q, Q and ([1]). It contains intractable 
quantities but can be accurately approximated given the time discretisation 
of the diffusion path. An alternative reparame trisation, defined in discrete 
time, was suggested in IChib and otherd ( 20061 ). It uses the transformation 



below 

xt - (xt_5 - 6fixixt,9x)) ^^ ,_. 

wt = — ,Vt. (5) 

ax{xt-s,Ox) 

In our setting the driving Brownian motion of Xq*^, denoted by Wq^^ and 

provided by ([5]), can be used to provide a discrete skeleton of the state 
vector 

Vo^^ = hUwtZxo,Ox,e,). (6) 

The model is now given by Wq'^^, that can be transformed to xt for which the 
Euler-Maruyama approximation can be used, and [T] which can be approxi- 
mated using the discretised state vector Vi^l^. 

Data augmentations schemes can be used for the models above. Gibbs 
versions of such schemes will alternate between updating Uq!* (or Wq*^) 
and consequently Vq^^* given 6, and 6 conditional on either Uq** (or Wq^^). 
The first step c a n be done either by the overlap ping block strategies in 



Chib and otherd (J2006l ) and iKalogeropoulosI (J2007l ) or with a particle filter 



in the context of a particle Gibbs algorithm. The second step of updating dx 
given Uq*^ or Wq^^ is usually implemented through a random walk Metropolis 
algorithm: 

• Let 6^ and Vi^l^'^ be the current values of 6x and Vjj:^^ respectively. 
Propose 6* from q(9*\9^). 

• Compute Vofjf* = huiuo;n,xo,9*,e^) 

• Accept with probability 

7r{9l,V,t'*\yi:n,ui^,^*9,,9y)c 



1 A 



A0'x,v,T\yi:n,<A'ev,9y)q{9i\ 



REFERENCES 



For the IChib and othersl ( 20061 ) formulation, Ug*^ can simply be replaced 



with Wq^^ in the algorithm above. 

Unfortunately both of the above algorithms may perform poorly. The 
problem is that every proposed value of 0^ implies a proposed trajectory 
of the ODE states vector Vt- As parts or functionals of this trajectory are 
observed with error, the proposed value of 0^ will not be accepted unless its 
associated Vt trajectory is close to these observations. Consequently, only 
small steps can be made on the 6^ space and the algorithm mixes very slowly. 
The problem intensifies as the noise variance becomes smaller and as the time 
horizon of the epidemic increases. Implementations of such algorithms in the 
simulated and real data of this paper are in line with this argument. Figure [6] 
displays the posterior draws for a in the dataset of the simulation experiment 
1 of Section 4 in the main paper. The posterior draws of a were obtained 
from a particle Gibbs algorithm combined with the algorithm above; note 
that in this model ujj*^ and it^g*^ are equal. In order to isolate the problem, 
the algorithm was run on a and /Jj only, and all the other parameters were 
held fixed at the values they where simulated from (a value of r = 0.1 was 
used). The 'true value' of a was 0.07 and we used a, 6 = 0.1. As clearly 
shown in the traceplot the mixing of the chain is quite poor, thus casting 
doubts on the reliability of its output. The difference in mixing quality with 
the corresponding traceplot of Fig. H] (bottom middle plot), corresponding 
to the PMCMC algorithm, is substantial. 

To sum up, both formulations of Gibbs data augmentation schemes (with 
or without reparametrisation) are very likely to lead to inaccurate and in- 
efficient MCMC algorithms. The use of particle MMH algorithms, termed 
as PMCMC in this paper, is therefore essential and the main paper focused 
on its implementation on diffusion driven epidemic models. PMCMC seems 
to provide a solution to the problem, but future research on Gibbs schemes 
with alternative reparametrisations addressing this problem, would be very 
helpful. 
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MSE 


Bias 


EKF 


0.0269 


0.0285 


Particle filter 


0.0242 


0.0064 


Improvement with regards to EKF 


-10% 


-77% 


Particle smoother 


0.0032 


0.0027 


Improvement with regards to P. filter 


-87% 


-64% 



Table 1: Mean Squarred Error and Bias of /3j estimates provided by the 
EKF, particle filter and particle smoother 





Exp l.a 


Exp l.b 


Exp l.b 


Exp 2.b 


r Simulation value 


0.1 


0.05 


0.1 


0.05 


Posterior mean 


0.103 


0.083 


0.078 


0.050 


Posterior median 


0.103 


0.084 


0.077 


0.050 


Posterior 95% c.i. 


[0.051; 0.152] 


[0.027; 0.137] 


[0.063; 0.96] 


[0.042; 0.060] 


a Simulation value 


0.07 


0.07 


n.d. 


n.d. 


Posterior mean 


0.066 


0.083 


0.016 


0.014 


Posterior median 


0.064 


0.084 


0.015 


0.014 


Posterior 95% c.i. 


0.048; 0.090 


0.046; 0.089 


0.010; 0.027 


0.001; 0.021 



Table 2: Mean, median and 95% confidence intervals for r and a estimates 
in four experiments 
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Posterior mean Posterior median 



Posterior 95% 



c.i. 



k-' 


1.56 


1.55 




[1.53 


1.60] 


t' 


1.00 


1.00 




0.92 


1.08 


/?cc(0) 


1.44 


1.36 




0.89 


2.30 


/3aa(0) 


1.40 


1.39 




1.25 


1.64 


/3ca 


0.30 


0.31 




0.16 


0.50 


Pac 


0.32 


0.32 




[0.18 


0.48] 


Ec{0) 


2.1 X 10"^ 


1.9 X 10-5 


[1.3 


X 10-5 


3.8 X 10-5] 


W) 


1.2 X 10^^ 


1.4 X 10-5 


[0.6 


X 10-5 


2.9 X 10-5] 


Rc{0) 


0.13 


0.12 




[0.06 


0.26] 


Ea{0) 


2.1 X 10-5 


2.0 X 10-5 


[1.1 


X 10-5 


3.5 X 10-5] 


4(0) 


1.0 X 10-5 


1.0 X 10-5 


[0.3 


X 10-5 


1.5 X 10-5] 


RaiO) 


0.16 X 10-5 


0.16 




[0.09 


0.28] 


0-c 


0.11 


0.10 




[0.08 


0.15] 


0-a 


0.08 


0.08 




[0.05 


0.11] 



Table 3: Mean, median and 95% confidence intervals for the parameters of 
the structured model applied to the A/HINI pandemic data 





r 


fc-l 


7-^ 


/3o 


E{0) 


m 


i?(0) 


0" 


2.5% quantile 


0.04 


1.55 


0.93 


0.80 


5.2 X 10-^ 


1.6 X 10-^ 


0.02 


0.04 


25% quantile 


0.09 


1.57 


1.03 


1.16 


1.6 X 10-5 


7.0 X 10-^ 


0.12 


0.05 


Median estimate 


0.11 


1.59 


1.08 


1.35 


2.3 X 10-5 


1.6 X 10-5 


0.17 


0.06 


75% quantile 


0.13 


1.60 


1.13 


1.56 


3.1 X 10-5 


2.8 X 10-5 


0.22 


0.07 


97.5% quantile 


0.17 


1.63 


1.23 


2.13 


5.2 X 10-5 


6.5x10-5 


0.33 


0.10 


Median when R(0) shifted +10% 


0.11 


1.59 


1.09 


1.41 


2.0 X 10-5 


1.9 X 10-5 


0.19 


0.06 


Median when R(0) shifted +20% 


0.11 


1.59 


1.09 


1.44 


1.8x10-5 


2.1 X 10-5 


0.24 


0.06 


Median when R(0) shifted -10% 


0.11 


1.59 


1.08 


1.31 


2.2 X 10-5 


2.2 X 10-5 


0.15 


0.07 


Median when R(0) shifted -20% 


0.12 


1.59 


1.09 


1.27 


1.9 X 10-5 


2.1 X 10-5 


0.13 


0.07 


Median when fc-^ shifted +10% 


0.11 


1.59 


1.09 


1.33 


1.8 X 10-5 


2.2 X 10-5 


0.15 


0.06 


Median when A;"^ shifted +20% 


0.12 


1.60 


1.08 


1.34 


1.8x10-5 


2.3 X 10-5 


0.17 


0.06 


Median when k~^ shifted -10% 


0.12 


1.58 


1.08 


1.28 


2.0 X 10-5 


2.0 X 10-5 


0.16 


0.06 


Median when k'^ shifted -20% 


0.12 


1.57 


1.09 


1.31 


1.9 X 10-5 


2.0 X 10-5 


0.16 


0.06 


Median when 7-^ shifted +10% 


0.11 


1.59 


1.12 


1.23 


2.2 X 10-5 


2.1 X 10-5 


0.15 


0.07 


Median when 7-^ shifted +20% 


0.10 


1.59 


1.14 


1.18 


2.0 X 10-5 


2.4 X 10-5 


0.15 


0.07 


Median when 7"-'^ shifted -10% 


0.11 


1.59 


1.06 


1.37 


1.9 X 10-5 


2.0 X 10-5 


0.16 


0.07 


Median when 7"^ shifted -20% 


0.10 


1.59 


1.02 


1.46 


2.0 X 10-5 


1.8 X 10-5 


0.17 


0.07 



Table 4: Original estimates compared to the ones resulting from respectively 
tilting the priors on R{0), 7-^ or k'^ by +10, +20, -10 or -20% 
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E[p{T|y)] (when t=0.05) 




3.5 1.3 0.47 0.17 0.064 

euler discr. timestep in days (logscale) 
E[p{T|y)] (when t=0.05) 




3.5 1.3 0.47 0.17 0.064 

euler discr. timestep in days (logscale) 



0.14 
0.12 



E[p(o|y}] (when t=0.05) 



3.5 1.3 0.47 0.17 0.064 

euler discr. timestep in days (logscale) 
E[p(G|y)] (when t=0.05) 



3.5 1.3 0.47 0.17 0.064 

euler discr. timestep in days (logscale) 



Figure 1: Convergence of the posterior density as the Euler discretization 
time-step 6 decreases (x-axis in the log-scale) 
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Figure 2: Acceptance rate as a function of N parts ^ in two situations where 
the noise amplitude is respectively 10% (full line) and 5% (dotted line). 
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Total Influenza Incidence 



Total Influenza Incidence 
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Figure 3: Illustration of how the underlying dynamic of the effective contact 
rate can be estimated from weekly recorded cases (r = 0.05) 
Green dots indicate simulated observed incidence (top panels). Green lines 
indicate simulated effective contact rate trajectories (bottom panels). Black 
dotted lines indicate the mean of the pointwise posterior density. Dark and 
light blue show credible intervals, respectively at 50% and 95% levels. 
Top left: experiment l.b, weekly number of cases observed with noise 
Top right: experiment 2.b, weekly number of cases observed with noise 
Bottom left: experiment l.b, simulated and estimated trajectory of the ef- 
fective contact rate 

Bottom right: experiment 2.b, simulated and estimated trajectory of the 
effective contact rate 
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Figure 4: MCMC traceplots for each component of 6 
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IBM "real-time" estimates 
DIG = 73.6 



r 



3 


20 Jul 


14 Sep 
DIC = 167 


9 Nov 
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Figure 5: Modeling choices and implications, aiming for robustness 
Black dotted lines indicate the mean of the pointwise posterior density. Dark 
and light blue show credible intervals, respectively at 50% and 95% levels. 
Left panels: estimates from an alternative modeling approach: exploring the 
full posterior density of an IBM diffusion model (left) 

Right panels: estimates from an alternative methodological approach: ex- 
ploring the posterior density of a BM diffusion model conditionned on a 
likelihood maximizing parameter 6* provided by the MIF algorithm (right) 
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Figure 6: MCMC traceplot for a when using a Particle Gibbs scheme with 
reparametrization 



